Effect of boundaries on the force distributions in granular media 
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The effect of boundaries on the force distributions in granular media is illustrated by simulations 
of 2D packings of frictionless, Hertzian spheres. To elucidate discrepancies between experimental 
observations and theoretical predictions, we distinguish between the weight distribution V(w) mea- 
sured in experiments and analyzed in the q- model, and the distribution of interparticle forces P(f). 
The latter one is robust, while V(w) can be obtained once the local packing geometry and P(f) are 
known. By manipulating the (boundary) geometry, we show that V[w) can be varied drastically. 
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A crucial property of granular materials is their he- 
terogeneity |Q. In particular, the strong fluctuations of 
interparticle forces and the organization of the largest of 
these in tenuous force networks have recently attracted 
considerable attention g |, |, |, §, g, g] . The probability 
density function of forces is thus a basic object of study. 
Measurements ^, 0, 0, p), numerical simulations || |7j 
and theory Q agree that such force distributions decay 
exponentially for large forces 0. The behavior for small 
forces is less well settled; while the g-model [|| seems to 
predict a vanishing probability, experiments and numeri- 
cal simulations clearly show that this probability remains 
non-zero for small forces. Since the small force distribu- 
tion may be a fingerprint of arching |^p|| , or of whether 
a system is jammed or unjammed 11 , it is important to 
obtain a clear physical interpretation of this discrepancy. 

In this paper, we resolve this issue by elucidating the 
effect of the local packing geometry on the force network 
in 2D packings of frictionless, Hertzian spheres under 
gravity (see Fig. la). To do so, it is important to dis- 
tinguish the effective weight W of a particle j from the 
interparticle forces F (Fig. lb). In the bulk, this weight 
is carried by the other particles on which particle j rests; 
however, the particles in the bottom layer are typically 
supported by the bottom only, so the forces exerted on 
the support are then to a good approximation the same 
as the weights. Thus it is essentially the distribution of 
weights V(W) which is probed in experiments where the 
particle-wall forces are extracted from the imprints on 
carbon paper [Q, || or by force sensors Q . Likewise, the 
main prediction of the g-model is for the distribution 
V{W)i rather than for the distribution of interparticle 
forces, denoted by P{F). 

For our case of frictionless spheres, we define the 
weights W as (see Fig. lb) 
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Here rrij denotes mass, g denotes gravity, Fy- are the 
interparticle forces, and the sum runs over all n c particles 




mg I 



(b) 




0.8 
(w) 
0.4 

0.2 
0.0 



•(c) , 


» ■ • s! 
* • • §i 


1.5 

<w 2 > 












• | 


1 1.0 










% 15 h 45 















1 



w 



FIG. 1: (a) Detail of a typical packing and force network in 
our simulations; the height h denotes the distance from the 
bottom, (b) Definition of interparticle forces F and weight 
W , for a particle with n c = 2. (c) The weight distribution 
V(w) at various heights between 10 and 30 in the bulk (open 
circles), for 2 < h < 3 (full curve) and at the bottom (dots). 
Inset: The second moment (w 2 ) as a function of height h. 



that exert a force on particle j from above (see Fig. lb) . 
In the following we rescale W and F to their average 
values, and write the rescaled weights and forces as w 
and /, with distributions V(w) and P{f). 

Our main findings are the following: (i) V(w) changes 
qualitatively when approaching a boundary; in particular 
the probability of finding a small weight is much larger at 
the bottom than in the bulk (Fig. lc). (ii) The number of 
contacts from above, n c , crucially influences the distribu- 
tion V(w), as can be anticipated from Eq. ([!]); the differ- 
ence between bulk and bottom "P(uj)'s is almost entirely 
due to the change in n c caused by the change in packing 
near the boundary (Fig. 2). (in) The force probability 



distributions P(f) and P(f z ) show a much weaker varia- 
tion when approaching the boundaries (Fig. 3). (iv) The 
distribution of n c 's near the bottom can be manipulated 
by, e.g., curving the boundary of a highly monodisperse 
packing and this can have a large effect on T{w) (Fig. 4). 

Numerical method Our 2D packings consist of fric- 
tionless particles under gravity; the particles interact 
through normal Hertzian forces, where / oc d 3 ^ 2 and d 
denotes the overlap distance Q. Unless noted other- 
wise, the material constants and gravity are chosen such 
that a particle deforms 0.1% under its own weight, and 
the particle radii are drawn from a flat distribution be- 
tween 0.4 < r < 0.6. Masses are proportional to the 
radii cubed. The container has a width of 24, employs 
periodic boundary conditions in the horizontal direction 
and has a bottom consisting of a fixed hard support. The 
data shown in this paper were obtained from 1100 real- 
izations containing 1180 particles each. We construct our 
stationary packings by letting the particles relax from a 
gas-like state by introducing a dissipative force that acts 
whenever the overlap distance d is nonzero. 

Distribution of weights In Fig. lc we show the weight 
distribution V(w) for the bottom particles (dots) which 
differs qualitatively from the bulk distributions (open cir- 
cles). Moreover, we observe that the transition is remark- 
ably sharp: in the slice 2 < h < 3, the weight distribution 
is already bulk-like (full curve). In the inset of Fig. lc we 
plot (w 2 ) which quantifies the width of V{w) as a func- 
tion of height h. The sharp transition of V(w) near the 
bottom is clearly visible. Additionally, in the bulk, (w 2 ) 
slowly increases with height, due to the deformations of 
the particles; this effect disappears for harder particles 
|ji~3"|| . Finally, near the top layer V(w) becomes sharply 
peaked and (w 2 ) decreases. 

To understand the change of V{w) near the bottom, 
consider the typical packing of Fig. 2a. The support 
aligns the bottom row of particles and thus strongly af- 
fects the directions of the interparticle forces. The forces 
between neighboring bottom particles are almost purely 
horizontal, so these hardly contribute to either the force 
on the particle from above [the "weight" in ([!])] or to 
the force needed to support it. This approximation be- 
comes better the smaller the polydispersity is. Thus, the 
average value of n c , the number of particles that press 
on the particle from above and hence contribute to its 
weight, is on average lower at the bottom than in the 
bulk. Intuitively it is clear that the probability of finding 
a small value of w increases with smaller n c for non- 
tensile forces. This statement can be made precise by 
considering Eq. (|l|) for fixed n c and analyzing V nc (w), 
the weight probabilities restricted to particles of given 
n c . As long as the joint probability distribution of the 
interparticle forces remains finite for small forces, it fol- 
lows from a phase-space argument that 
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FIG. 2: (a) Detail of a typical packing, showing the domi- 
nance, near the bottom, of layer-to- layer forces (black lines) to 
intralayer forces (white lines) in determining w. The numbers 
show the values of n c for the respective bottom particles. (b,c) 
Decomposition of V{w) according to Eq. (^) in the bulk (b) 
and at the bottom (c); The measured bulk values for the frac- 
tions {po, pi, P2, P3} in Eq. (3) are {0.01, 0.11, 0.52, 0.36}, and 
the bottom values are {0.08,0.46,0.44,0.02}; at the bottom, 
we excluded the intralayer, almost horizontal forces. (d,e) 
When rescaled to the average value for each distribution func- 
tion, Vi(w) (d) and V^w) (e) are essentially the same in the 
bulk (open circles) and at the bottom (dots). 



for all n c > 1. The particles which do not feel a force 
from above, n c = 0, give a 5— like contribution at W = 1; 
for deep layers this occurs for w<l. 

To check the validity of this idea, we have determined 
P na (w) both in the bulk and near the bottom by deter- 
mining n c for each particle and decomposing the weight 
distribution P(w) into the ■p nc (w;)'s, 
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The crucial difference between bottom and bulk are the 
fractions p Uc of particles that feel n c other particles press- 
ing on them from above. This can be seen in Figs. ||b- 
c, where we show the decomposition of V(w) according 
to Eq. (|J). This picture is also confirmed by Figs. ||d- 
e, which show that the individual distribution functions 
Pi(w) and V2 (w) do not differ significantly between bulk 
and bottom (to compare these, we have to normalize 
them not to the total average weight in each layer, but 
to the average weight of the particles with the same n c ). 
Note also that Figs. §b-e show that Eq. (Q) is valid, ex- 
cept for V\{w) at the bottom for small values of w; this 
deviation comes from neglecting the intralayer, "almost 
horizontal" forces, whose small vertical components even- 
tually make P\(w) — > for very small weights. 
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FIG. 3: (a) P(f) in the bulk (open circles) and for the layer- 
to-layer forces (see Fig. 2a) near the bottom (dots), (b) P(f z ) 
in the bulk (open circles) and for the layer-to-layer forces near 
the bottom (dots). The solid line is obtained by integrating 
P(f) over all angles (see text). (c,d) Scatter plot of (/y, <fij) 
for (c) the bulk forces and (d) the layer-to-layer forces near 
the bottom. Same packings as for Figs. |l] andg. 

Interparticle forces The distributions of the interpar- 
ticle forces are much more robust than the weight dis- 
tributions. The results for the distribution P(f) of |/| 
are shown in Fig. [}]a. The only difference between bulk 
and bottom distributions is that the small peak around 
/ = 0.7 for bulk forces becomes a plateau for P(f) for 
the forces close to the bottom (see Fig. 2a). It is intrigu- 
ing to note that this change from a plateau to a peak 
is reminiscent of what is proposed as an identification of 
the jamming transition . 

We also have measured the distribution of angles <py, 
which define the orientation of the , and find that these 
angles are uniformly distributed and independent of the 
absolute value of / in the bulk, see Fig. ||c. Thus, in the 
bulk our packing is isotropic. Near the boundary, how- 
ever, this isotropy is broken strongly: in agreement with 
our scenario for the influence of the packing geometry, the 
angles of the forces between bottom particles and those in 
the layer above are concentrated around ir/3 and 2ir/3, 
as Figs. ||a and ||d show. Near the bottom, therefore, 
the interparticle forces naturally divide up into almost 
horizontal intralayer forces and "layer-to-layer" forces. 

Since the weight distribution is determined by the z- 
components of the forces, let us also investigate the dis- 
tribution P(f z ). According to Fig. ||b, P(f z ) also re- 
mains non-zero for small f z , both in the bulk and near 
the bottom. There is a substantial difference, however, 
associated with the difference in packing. In the bulk, we 
have seen in Fig. ||c that there is no noticeable correlation 
between the angles ipy and the force strength /y . Hence 
in the bulk P(f z ) « Jd<pdf P(f)P(ip) S(f z - fain(<p)) 
with P(ip) = const. Indeed, the distribution obtained 
by numerical integration of this relation with P(f) from 
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FIG. 4: (a,b) Packing and force networks in a weakly poly- 
disperse packing near a flat bottom (a) and a curved bot- 
tom (b). (c) Distribution of weights V(w) on the flat bottom 
(solid line), convex curved bottom (dashed line) and the con- 
cave curved bottom (dotted line). The various shapes orig- 
inate from the corresponding {po,pi,p2}: {0.00,0.10,0.90}, 
{0.02, 0.39, 0.58} and {0.04, 0.46, 0.50} respectively. 

Fig. |^a and a uniform angle distribution yields the solid 
line in Fig. ||b, which closely follows P(f z ) as measured in 
the bulk. Near the bottom, on the other hand, the value 
of sin(ip) is concentrated around \V2> » 0.866: in the 
approximation that the distribution of sin(y) is sharply 
peaked at this value, the shape of P(f z ) is close to that 
of P(f), which is confirmed by direct comparison of the 
dotted datasets of panels a and b of Fig. |[ 

Manipulating V(w) In the previous paragraphs we 
have shown that the weight distribution V(w) is very 
sensitive to the local packing geometry, while the dis- 
tribution of interparticle forces is robust. This allows 
one to manipulate V(w) at the bottom by changing the 
boundary conditions. We illustrate this by simulations 
of weakly polydisperse particles, 0.49 < r < 0.51, both 
with a flat bottom (Fig. |^a) and a curved bottom, con- 
sisting of two circle segments of radius 20 glued together 
(Fig. 0b). For the flat bottom, the particles form an al- 
most perfect hexagonal packing, leading to particles with 
mostly n c = 2 (horizontal contacts again not included). 
In agreement with (ph, V(w) increases linearly for small 
w in this case. The weakly curved bottom locally dis- 
turbs this crystalline structure (Fig. 4b), causing a dra- 
matic change in the fractions p„ c , and correspondingly in 
V(w) (Fig. |]c). Note that the small difference between 
the distributions of the convex and the concave part of 
the bottom are reflected in the p n as well. Interestingly, 
P(f) is in all these cases indistinguishable from P(f) in 
the strongly polydisperse case |l3| . 

Perspective The message that emerges naturally from 
the above analysis is clear: in experiments in which the 
forces on a boundary are probed, one measures effectively 
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the weights w. These weights, however, are not the most 
fundamental quantities of a granular packing, as they are 
derived from the interparticle forces fy. These capture 
the full microscopic structure, and the distribution func- 
tion P(f) of the force strength is quite insensitive to the 
packing, in contrast to V(w). Simple phase space consid- 
erations show that a grain with n c contacts with particles 
that press downwards on it makes a contribution to the 
weight distribution V(w) which scales as w n "~ 1 as w — > 0. 
Thus, the small weight distribution is dominated by par- 
ticles with few such contacts, in particular by those with 
just one, and the change in geometry near boundaries 
leads to an atypical V(w). In addition, we occasionally 
observe a small peak at w <C 1 due to "loose" particles 
(n c — 0) which do not feel a force from above; there are 
indications for such a peak of V(w) at w = in recent 
precise experiments (see Fig. 5 of ||). 

In the standard g-model, weights are randomly dis- 
tributed over a fixed number of neighbors one layer be- 
low; in the simplest version there are two such neighbors 
that receive a fraction q and 1 — q of the weight || . Due 
to a fixed connectivity, this model cannot be expected 
to capture the behavior of V(w), especially near bound- 
aries. The product of the weight w and q, however, could 
be interpreted as an interparticle force. Interestingly, in 
the simplest case of a uniform probability distribution 
of the q's, the probability distribution P(qw) is a pure 
exponential (ll| |l4) . In simulations of the q- model with 
random connectivities, a variety of "P(u>)'s can be ob- 
tained, similar to what we found here for the frictionless 
spheres p3[ . 

It is possible to test our framework in experiments. For 
P(w) we expect that curvature effects, such as shown in 
Fig. 4, only play a role when they break highly ordered 
packings; indeed carbon-paper measurements at the side- 
walls in cylinders give a very similar V(w) as near the 
bottom H. A simpler way to change the boundary con- 
ditions may be to include a layer of larger particles at 
the boundary; their value of n c will be higher, and we 
expect that V(w) for small w will decrease for larger n c . 
Furthermore, there is a strong need of direct determina- 
tion of P(f), both in the bulk and near the boundaries, 
since the present data concerns V{w) only |3[ |j. An 
interesting observation in the context of "jamming" fTlfl 
is that in our simulations the peak in P(f) appears to 
vanish near the boundaries (Fig. 3a); are granular mate- 
rials no longer jammed here, and is this relevant for the 
localization of shear bands near boundaries? 

An important issue for future study is clearly the role 
of friction and dimensionality. Our numerical study has 
been done in two dimensions with frictionless spheres; 
however, recent studies indicate M that the coordination 
number for 3D packings with friction is similar to those 



have advanced is therefore expected to capture the real- 
istic case of three dimensions with friction, because our 
phase space arguments are independent of dimension. 

We finally note an interesting open issue. In the q- 
model, V(w) approaches its asymptotic expression as a 
function of depth algebraically slow |L4j). In our pack- 
ings, the convergence appears to be much faster — is this 
a real discrepancy and if so is it related with the same 
packing issues? 

We are grateful to Martin Howard and Hans van 
Leeuwen for numerous illuminating discussions. 
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